trait <- read_csv("data/dummy_trait.csv")Rows: 77 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): sp
dbl (8): lma, ll, a_mass, r_mass, n_mass, p_mass, wd, sm
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.tree <- read.tree("data/dummy_tree.newick")x <- seq(1, 6)
y <- seq(1, 10)
z <- matrix(numeric(length(x) * length(y)), ncol = length(y))
z[1,8:9] <- z[2,8:9] <- z[3,8:9] <- 1
z[4,2:4] <- z[5,2:4] <- z[6,2:4] <- 1
z2 <- matrix(numeric(60), ncol = length(y))
z2[,1:5] <- 1
env <- tibble(env1 = as.numeric(z),
env2 = as.numeric(z2)) %>%
mutate(row = rep(1:6, 10)) %>%
mutate(col = rep(1:10, each = 6)) %>%
mutate(site = str_c("site", 1:60)) %>%
mutate(topo = ifelse(env1 == 1, "Ridge", "valley")) %>%
mutate(type = ifelse(env2 == 1, "Rainforest", "Rubber"))
ggplot(env, aes(y = row - 0.5, x = col - 0.5, fill = topo)) +
geom_raster(alpha = 0.8) +
geom_vline(xintercept = 0:10) +
geom_hline(yintercept = 0:6)
ggplot(env, aes(y = row - 0.5, x = col - 0.5, fill = type)) +
geom_raster(alpha = 0.8) +
geom_vline(xintercept = 0:10) +
geom_hline(yintercept = 0:6)
n_sp <- nrow(trait)
#beta_LMA <- rnorm(n_sp, log(LMA), 0.2)
beta_LMA <- 1.5
beta_WD <- 3
ts <- trait[,-1] %>% log %>% scale %>% as_tibble11 * 7
generate_env <- function(seed = 123, model = c("ridge_negative", "ridge_positive", "random")) {
set.seed(seed)
x <- seq(1, 13)
y <- seq(1, 9)
z <- matrix(rnorm(length(x) * length(y)), ncol = length(y))
if (model == "ridge_positive") {
z[3:6, 6:9] <- 2
z[9:11, 1:4] <- 2
} else if (model == "ridge_negative") {
z[3:6, 6:9] <- -5
z[9:11, 1:4] <- -5
}
for (i in 2:(length(x)-1)) {
for (j in 2:(length(y)-1)) {
z[i, j] <-
rnorm(1,
mean(
c(z[i - 1, j - 1],
z[i - 1, j],
z[i - 1, j + 1],
z[i, j - 1],
z[i, j + 1],
z[i + 1, j - 1],
z[i + 1, j],
z[i + 1, j + 1])),
0)
}
}
z <- z[,-9]
z <- z[,-1]
z <- z[-13,]
z <- z[-1,]
x <- seq(1, 11)
y <- seq(1, 7)
x2 <- rep(x, length(y))
y2 <- rep(y, each = length(x))
# as.numeric(z)
# as.numeric(t(z))
dat <- tibble(x = x2 |> as.numeric(),
y = y2,
z = as.numeric(z) )
# mutate(z = ifelse(x > 5, z + 1, z))
dat2 <- dat %>%
mutate(x = x - 1) %>%
mutate(y = y - 1)
dat2
}
moist <- generate_env(model = "ridge_negative")
nitro <- generate_env(model = "random")
phos <- generate_env(model = "random")
nitro <- nitro |>
mutate(soil_n = exp(z)) |>
mutate(soil_n = soil_n / max(soil_n) * 0.4)
phos <- phos |>
mutate(soil_p = exp(z)) |>
mutate(soil_p = soil_p / max(soil_p) * 0.2)
moist <- moist |>
mutate(soil_moist = z)
dat2 <- nitro |>
dplyr::select(-z) |>
mutate(soil_p = phos$soil_p) |>
mutate(soil_moist = moist$soil_moist)
ggplot(dat2, aes(x = x, y = y, fill = soil_n |> log())) +
geom_raster()
ggplot(dat2, aes(x = x, y = y, fill = soil_p |> log())) +
geom_raster()
ggplot(dat2, aes(x = x, y = y, fill = soil_moist)) +
geom_raster()
tmp_fun <- function(data, env) {
tmp <- NULL
for (i in 0:9) {
for (j in 0:5) {
s1 <- data |>
filter(x == i & y == j) |>
pull({{env}})
s2 <- data |>
filter(x == (i + 1) & y == j) |>
pull({{env}})
s3 <- data |>
filter(x == i & y == (j + 1)) |>
pull({{env}})
s4 <- data |>
filter(x == i & y == (j + 1)) |>
pull({{env}})
tmp <- c(tmp, mean(c(s1, s2, s3, s4)))
}
}
tmp
}
env3 <- env |>
mutate(soil_moist = tmp_fun(dat2, soil_moist))
ggplot(env3, aes(x = col, y = row, fill = soil_moist)) +
geom_raster()
dat2 |>
dplyr::select(x, y, soil_n, soil_p, soil_moist) |>
round(3) |>
write_csv("data/soil.csv")Abundance of ith species at site j:
\[ Abundance_{ij} \sim Pois(\lambda_{ij}) \] \[ \log\lambda_{ij} = \beta_{0j} + \beta_{1j} \times Ridge_i + \beta_{2j} \times Forest_i + \beta_{3j} \times Ridge_i \times Forest_i \] \[ \beta_{mj} = \gamma_{m0} + \gamma_{m1} \times logLMA_j + \gamma_{m2} \times logWD_j + \epsilon_{mj} \] \[ \epsilon_{mj} \sim N(0, \sigma_m) \]
beta0 <- rnorm(n_sp, 5, 0.2)
beta1 <- rnorm(n_sp,
mean = -5 + 1.5 * ts$lma + 0.6 * ts$wd,
sd = 0.2)
beta2 <- rnorm(n_sp, -1 + 0.8 * ts$sm, 0.2)
beta3 <- rnorm(n_sp,
mean = 1 * ts$lma,
sd = 0.2)
beta4 <- rnorm(n_sp,
mean = -0.5 * ts$lma,
sd = 0.2)
beta <- cbind(beta0, beta1, beta2, beta3, beta4)
env2 <- tibble(intercept = 1,
env1 = env$env1,
env2 = env$env2,
interaction = env$env2 * env$env1,
moist = env3$soil_moist)
env2# A tibble: 60 × 5
intercept env1 env2 interaction moist
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0 1 0 0.193
2 1 0 1 0 0.0693
3 1 0 1 0 -0.280
4 1 0 1 0 -1.18
5 1 0 1 0 -1.86
6 1 0 1 0 -2.41
7 1 0 1 0 0.198
8 1 0 1 0 0.0913
9 1 0 1 0 -0.652
10 1 1 1 1 -1.94
# ℹ 50 more rowssummary(ts) lma ll a_mass r_mass
Min. :-2.35492 Min. :-1.9081 Min. :-2.63320 Min. :-2.86039
1st Qu.:-0.77439 1st Qu.:-0.6445 1st Qu.:-0.51068 1st Qu.:-0.65759
Median :-0.07205 Median :-0.1391 Median :-0.01897 Median : 0.03553
Mean : 0.00000 Mean : 0.0000 Mean : 0.00000 Mean : 0.00000
3rd Qu.: 0.73571 3rd Qu.: 0.4938 3rd Qu.: 0.64283 3rd Qu.: 0.64807
Max. : 2.91324 Max. : 2.8215 Max. : 2.39729 Max. : 2.05776
n_mass p_mass wd sm
Min. :-2.7071 Min. :-2.65152 Min. :-3.20202 Min. :-2.14566
1st Qu.:-0.8276 1st Qu.:-0.69644 1st Qu.:-0.73935 1st Qu.:-0.65216
Median : 0.1815 Median : 0.05989 Median :-0.01349 Median :-0.02183
Mean : 0.0000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
3rd Qu.: 0.7325 3rd Qu.: 0.67051 3rd Qu.: 0.62138 3rd Qu.: 0.58372
Max. : 2.1057 Max. : 1.92002 Max. : 2.15467 Max. : 2.51892 mu <- beta %*% t(env2)
ab <- matrix(numeric(60*77), ncol = ncol(mu))
for (i in 1:nrow(mu)) {
for (j in 1:ncol(mu)) {
ab[i, j] <- rpois(1, exp(mu[i, j]))
#ab[i, j] <- rpoilog(1, exp(mu[i, j]), 1)
}
}# Function to plot abundance for top/bottom species based on a specific trait
plot_abundance <- function(env, ab, species_list, title_list) {
for (i in seq_along(species_list)) {
p <- env %>%
mutate(abund = ab[species_list[i], ]) %>%
ggplot(aes(x = col - 0.5, y = row - 0.5, fill = abund)) +
ggtitle(title_list[i]) +
geom_raster() +
geom_vline(xintercept = 0:10) +
geom_hline(yintercept = 0:6)
print(p)
}
}
# Set rownames for abundance matrix `ab`
rownames(ab) <- trait$sptop_sp_lma <- trait %>%
arrange(-lma) %>%
slice(1:5) %>%
pull(sp)
plot_abundance(env, ab, top_sp_lma, top_sp_lma)
low_sp_lma <- trait %>%
arrange(lma) %>%
slice(1:5) %>%
pull(sp)
plot_abundance(env, ab, low_sp_lma, low_sp_lma)
top_sp_sm <- trait %>%
arrange(-sm) %>%
slice(1:5) %>%
pull(sp)
plot_abundance(env, ab, top_sp_sm, top_sp_sm)
low_sp_sm <- trait %>%
arrange(sm) %>%
slice(1:5) %>%
pull(sp)
plot_abundance(env, ab, low_sp_sm, low_sp_sm)
abundant_sp <- c(top_sp_lma, low_sp_lma, top_sp_sm, low_sp_sm) |>
unique() |>
sort()
abundant_sp [1] "Beilschmiedia_robusta" "Camellia_forrestii"
[3] "Camellia_Pitardii" "Castanopsis_wattii"
[5] "Claoxylon_khasianum" "Eurya_trichocarpa"
[7] "Ficus_henryi" "Ficus_heteromorpha"
[9] "Itoa_orientalis" "Lindera_caudata"
[11] "Lithocarpus_crassifolius" "Lithocarpus_henryi"
[13] "Lithocarpus_xylocarpus" "Machilus_salicina"
[15] "Neolitsea_polycarpa" "Rhododendron_decorum"
[17] "Rhus_sylvestris" "Schima_noronhae"
[19] "Symplocos_anomala" plot_abundance(env, ab, abundant_sp, abundant_sp)
rare_sp <- trait %>%
filter(!(sp %in% abundant_sp)) %>%
pull(sp) |>
# setdiff(trait$sp) |>
sort()
rare_sp [1] "Acer_campbellii" "Actinodaphne_forrestii"
[3] "Alnus_nepalensis" "Anneslea_fragrans"
[5] "Betula_alnoides" "Camellia_assamica"
[7] "Camellia_yunnanensis" "Castanopsis_echidnocarpa"
[9] "Craibiodendron_yunnanense" "Cyclobalanopsis_stewardiana"
[11] "Eriobotrya_bengalensis" "Eriobotrya_tengyuehensis"
[13] "Eurya_paratetragonoclada" "Eurya_quinquelocularis"
[15] "Exbucklandia_populnea" "Ficus_racemosa"
[17] "Glochidion_a" "Ilex_corallina"
[19] "Ilex_manneiensis" "Illicium_macranthum"
[21] "Laurocerasus_undulata" "Leucosceptrum_canum"
[23] "Lindera_thomsonii" "Lithocarpus_confinis"
[25] "Lithocarpus_dealbatus" "Lithocarpus_hancei"
[27] "Lithocarpus_truncatus" "Litsea_a"
[29] "Litsea_chunii" "Litsea_elongata"
[31] "Machilus_gamblei" "Machilus_yunnanensis"
[33] "Mahonia_duclouxiana" "Manglietia_insignis"
[35] "Melia_toosendan" "Meliosma_arnottiana"
[37] "Meliosma_kirkii" "Michelia_floribunda"
[39] "Myrsine_semiserrata" "Neolitsea_chuii"
[41] "Olea_rosea" "Osbeckia_mepalensis"
[43] "Padus_napaulensis" "Rhododendron_irroratum"
[45] "Rhododendron_leptothrium" "Saurauia_napaulensis"
[47] "Schefflera_fengii" "Schefflera_shweliensis"
[49] "Schima_argentea" "Skimmia_arborescens"
[51] "Sterculia_nobilis" "Stewartia_pteropetiolata"
[53] "Stewartia_yunnanensis" "Symplocos_grandis"
[55] "Symplocos_ramosissima" "Symplocos_sumuntia"
[57] "Ternstroemia_gymnanthera" "Vaccinium_duclouxii" tibble(sp = abundant_sp, abund = "abundant") |>
bind_rows(tibble(sp = rare_sp, abund = "rare")) |>
mutate(sp2 = str_replace_all(sp, "_", " ")) |>
write_csv("data/note.csv")plot_abundance(env, ab, rare_sp, rare_sp)
